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ABSTRACT 


The General Motors Research Laboratories has developed an 
image processing system which automatically analyzes the 
size distributions in fuel spray video images. Images are 
generated by using pulsed laser light to freeze droplet 
motion in the spray sample volume under study. This coherent 
illumination source produces images which contain droplet 
diffraction patterns representing the droplets degree of 
focus. Thousands of images are recorded per sample volume 
to get an ensemble average of the distribution at that spray 
location. After image acquisition the recorded video frames 
are replayed and analyzed under computer control. 

The analysis is performed by extracting feature data 
describing droplet diffraction patterns in the images. This 
allows the system to select droplets from image anomalies 
and measure only those droplets considered in focus. The 
system was designed to analyze sprays from a variety of 
environments. Currently these are an ambient spray chamber, 
a high pressure high temperature spray facility, and sprays 
in a running engine. 


Unique features of the system are the totally automated 
analysis and droplet feature measurement from the grayscale 
image. Also it can distinguish non- spherical anomalies from 
droplets which allows sizing of droplets near the spray 
nozzle. 
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This paper describes the feature extraction and image 
restoration algorithms used in the system. Preliminary 
performance data is also given for two experiments. One 
experiment gives a comparison between a synthesized 
distribution measured manually and automatically. The 
second experiment compares a real spray distribution 
measured using current methods against the automatic system. 
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INTRODUCTION 

Basic research is being conducted to relate the combustion 
process and the design and placement of fuel injectors. This 
is being accomplished by studying fuel spray droplet 
dynamics. The fuel spray study utilizes a system which 
records video images of spray droplets directly from a 
variety of fuel spray environments. These include an ambient 
spray chamber[l,2], a high temperature, high pressure test 
chamber, and the combustion chamber of a running engine[3]. 

Because of the harsh spray environments the analysis of 
video images at General Motors Research Laboratories is 
currently done manually by observing the images on a TV and 
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selecting only "in focus" droplets for velocity and size 
measurements. Manual data reduction is time consuming, 
tedious and inconsistent for different operators. These 
problems led to the development of an automatic system for 
droplet measurement. This system should access archived 
image data and extract droplet size information without 
manual intervention. The system must also analyze poor 
quality images containing low frequency intensity gradients 
which exceed droplet contrast levels. The realization of 
these requirements has pushed the capabilities of this spray 
analysis system beyond any other direct imaging 
system[4,5,6,7,S] . 


HARDWARE 

Figure 1 shows a block diagram of the hardware utilized in 
this study. For further details, refer to Oberdier 
paper [9]. Images are created by direct imaging of a spray 
sample area on a high resolution CCD (charge coupled device) 
array camera . This camera has a sensor resolution of 380 X 
480 picture elements (pixels). Illumination of the sample 
volume is provided by a pulsed 100 mJ nitrogen laser at a 
wavelength of 337 nm. The 10 ns laser pulses are collimated 
and passed through the test volume perpendicular to the 
spray. The imaging lense used has a 117 mm focal length. The 
recording camera is positioned to give a resolution of 
3.3 urn per image pixel. This magnification allows a 1.2 X 
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1.5 mm field of view at the focal plane. Because of the 
filter operation performed by the software to remove low 
frequency image degradation, the measurable droplet diameter 
range is limited at 5 to 50 pixels. This translates an 
effective droplet diameter range of 16 to 165 urn. Thus the 
imaging lense is switched to accommodate other ranges. 

Image acquisition is accomplished by using a microcomputer 
based controller which synchronizes laser firing, camera 
scanning and image storage to the experimental process. The 
controller was developed at General Motors and will 
synchronize image acquisition at a particular engine crank 
angle or free run at standard video rates (30 frames per 
second). Images are stored on a magnetic video disk recorder 
which allows full frame, random access of the images. After 
the data is collected the stored images are digitized and 
analyzed automatically by a Vicom image analysis system. 
This system digitizes each image to a 512 X 512 pixel array 
at 8 bits of resolution per pixel. It should be noted that 
current camera and recorder technologies limit the 
realization of this specification to 6 bits of of intensity 
resolution and the horizontal bandwidth to 450 lines. 
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IMAGE PROCESSING ALGORITHMS 


Image processing algorithms were applied to perform the 
following steps. 

1. Image normalization - Remove low frequency 

degradation in the recorded images. 

2. Segmentation - Find places in the image which may be 
the centroid of droplets. 

3. Feature extraction - Extract feature data from the 

image restored in step one at those centroid 

locations . 

4. Classification - Decide if the object is an in focus 
droplet by using extracted diffraction feature data. 

Image Restoration 

The fuel spray images recorded from the high temperature, 
high presure chamber and the running engine chamber have low 
frequency intensity degradation. This is caused by droplets 
which hit the observation window, refracting gradients in 
the optical path, and sensor scan variations. This 
degradation has intensity values which can be similar to 
droplet intensities. This degradation is removed in the 
following manner. 

As discussed in[10,ll,12], an image formation model can be 
represented mathematically as: 
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( 1 ) 


g(x,y) « d(x,y)f (x,y) ♦ n(x,y) ♦ b 


where 

g(x,y) * detected Image 
f<x,y) » original image 
d(x,y) = multiplicative noise 

u\a — auui Live uuxae 

b ■ intensity bias. 

We assume that the low frequency degrading function is represented 
by d(x,y) . The additive noise 

n(x,y) represents high frequency ( pixel to pixel ) 
digitization and camera sensor noise. The b term is added 
to the model by the authors to account for any intensity bias which 
may be added to the image by video circuitry or ambient light. 

Solving for the original image, f(x,y), before the degradation 
occurred results in: 

(2) f(x,y) » [g(x r y) - n(x,y) - b] / d(x,y). 


Instead of estimating the additive noise and subtracting it 
from the image g(x,y) , the image g(x,y) is passed through a 
low pass filter. This is accomplished by convolving g(x,y) 
with a Gaussian impulse response whose coefficients are 
given as 
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Figure 2 shows a 

plot 

of g(x,y) 

before 

and after 

the 

Gaussian operation. The intensity bias, b 

is calculated as 

being the lowest 

intensity value 

in the 

image after 

the 

additive noise is 

removed. This 

bias is 

subtracted 

from 

every image point 

in the 

image. 





The last step in this image normalization process is the 
removal of the low frequency degradation, d(x,y). This is 
first estimated by using a nonlinear filtering technique 
called morphological filtering. This was selected over 
other low pass filters because it preserves edge structure. 
The effect of this filter operation is to delete image 
objects smaller than the defined size of the filter. This 
leaves am image which contains only the background intensity 
levels. Therefore to estimate d(x,y), the filter size is 
selected to be larger than the largest object to be 
measured. 
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A simple representation of this filter as given in[13] is: 


(3) d(s,t) = max £ f(x,y) . circ ^(x-s)*+(y-t) 1 /r,33 

where d(s,t) = resulting image 
f(x,y) = original image 

circ [^/(x-s)*+(y-t) 2 /r 0 ] describes the circular 

structuring element with 
radius r 0 and equals: 

1 for >^x-s) 2 +(y-t) 2 /r 0 <= 1 

0 otherwise 

The max function propagates local image intensity maxima 

over the filter size defined by r, for every pixel in the 

image. The original boundaries are preserved by performing 
a reverse propagation on the transformed image: 

(4) d(s,t)= mint f(x,y) . circ [>/(x-s) 2 +(y-t) 2 /r,]] 

The min function propagates local image intensity minima 

over the filter size defined by r, for every pixel in the 

image. To normalize the image, the estimated d(x,y) is 
divided into the intermediate image, ( g(x,y) - n(x,y) - b) 
per equation (2) . 


Figure 3, 4 and 5 show an original image, its background 
estimate, d(x,y) using this filter, and the normalized image 
after division. 

SEGMENTATION - Locating candidate objects 

By utilizing the previous image restoration algorithm, 
intensity thresholding can be used to segment areas which 
may or may not be droplets, i.e., candidates [14] . Image 
threshold levels are selected by analyzing the intensity 
histogram [15, 16] . It is assumed that the largest 
distribution in the histogram corresponds to the background 
intensity. The midpoint of this distribution's positive 
slope is selected as the threshold level if the histogram is 
unimodal. If there is a peak prior to the background peak 
the valley between both distributions is selected. 


The output of the threshold operation is a binary image in 
which the dark areas correspond to candidates. Region 
boundaries are generated by using a 4 adjacency border 
algorithm [ 17 ] . This algorithm allows a fast way to generate 
all object boundaries at one time by treating the image as a 
whole. It is defined mathematically as 
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(5) B(x,y) = S(x,y) 


* [ S(x-l,y) + S(x+l,y) + S(x,y+1) + S(x,y-1) ] 


where 

B(x,y) * boundary image 

S(x,y) = original threshold image 

S(x,y) = logical negation of S(x,y) 

* = logical AND 

+ = logical OR. 

The output of this algorithm is object edge pixels which 
define the object's boundary. Figures 6 and 7 show a 
threshold and its associated boundary image. 

The image is scanned by a tracking algorithm which follows 
object boundaries and calculates the centroid, radius (a 
circle is fit to the points) and number that represents how 
good the boundary fits the circle. The centroids are points 
in the normalized image g(x,y) which will be used as feature 
extraction locations. 

Feature Extraction 

To extract features from candidate objects two curves are 
determined. Both curves are calculated by assuming symmetry 
around the centroids determined in the previous section. One 
curve is a radial intensity profile, the other a radial 
standard deviation. To accomplish this, the vector distance 
from each candidate's centroid to all neighboring points is 
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calculated. Neighboring points are averaged for an area 
spanning twice the radius value found when calculating the 
boundary centroid. The distance to neighboring points is 
modified to compensate for the aspect ratio induced by the 
camera. This is calculated using 0.8 times the vertical 
displacement. Intensity values for similar distance vectors 
are accumulated to get average radial intensity and radial 
standard deviation. See Figure 8. 

Because of the way the centroid was selected two things are 
possible. First, the calculated intensity and deviation 
curves may or may not represent a droplet. Many kinds of 
anomalies are possible. For example the object could 
represent a multiple droplet cluster, a blob that was 
created by the selected threshold level or an anomaly 
created by edge effects after creating an image mosaic. 
Examples of these are shown in Figure 7 . These anomalies 
are distinguished using the circular fit number generated 
during boundary tracking along with the intensity and 
deviation curves during the classification process. 
Secondly, the centroid may indicate a valid droplet but may 
not be the real object center. To correct this, the software 
attempts to reposition all droplet centers by using radial 
deviation information as an indication of best center 
position. Repositioning is done by calculating the standard 
deviation for the four quadrants of a candidate. See 
Figure 9. The derivative of the quadrant curves 2,3,4 is 
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cross correlated with the derivative of quadrant 1. The 
resulting correlation peaks are at index points which weight 
the direction of center repositioning in the following way: 

(6) Delta X = TRUNCATE ( 0.5 + cos* ( -Q2 - Q3 + Q4 )) (7) 
Delta Y = TRUNCATE ( 0.5 + sin* ( Q2 - Q3 - Q4 )) 

where * = 45 degrees =(.707) 

Ql-4 = correlation index values 


This procedure iterates up to 5 times per candidate before 
aborting if the delta valves do not converge to zero. 
Assuming that the radial intensity profile corresponds to 
the diffraction pattern of a real droplet, a variety of 
features describing the pattern are calculated. Some of 
these features are : 

1. Intensity profile slope at the droplet center. 

2. Number of rings within the droplet. 

3. Slope of intensity profile at droplet edge. 

4 . Droplet contrast . 

5. Intensity profile overshoot height and width. 

6. Droplet radius measured at the maximum profile slope 
at the droplet's edge. 

Extracted features and where the feature is measured is 
shown on the curves in Figure 10. This figure was generated 
using theoretical radial intensity plots of an ideal opaque 
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100 micron sphere at 200 and 400 microns from the plane of 
focus. 

A rudimentary classification procedure uses these and other 
features to select the radial intensity curves which 
represent droplets that are in the sample volume. 

Experimental Procedure 

Two preliminary experiments were run to quantify the 
accuracy of the automatic system vs . manual techniques . In 
the first, a slide of polystyrene microspheres was measured 
using a micrometer attached to a microscope (labeled 
MICROMETER in Figure 11). 500 measurements were made on the 
slide to get an ensemble average of the distribution. In a 
similar fashion 50 images at random slide positions were 
analyzed by the automatic system (VIDEO-AUTO) . Images were 
acquired by illuminating the slide with a 5 mw helium neon 
laser. A CCD camera inline with the laser detected a sample 
area using the same 'experimental' optics described in the 
hardware section. The slide was held in a micropositioning 
device and images with a high droplet density were recorded. 
These same images were also measured manually on a video 
monitor using a scale (VIDEO-MANUAL). 

In a second experiment, 400 images were recorded in the high 
pressure high temperature test facility using a Heptane 
gasoline spray. A nitrogen laser and vidicon camera was 
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used as the illumination source and sensor. The optics used 
were the same as in experiment 1. The images were manually 
measured on a video monitor by two different observers 
(VIDEO-MANUAL #1,#2) and also by the automatic system. 
Because this experiment measured a real spray the sample 
counts were corrected for depth of field effects. This 
procedure compensates for the fact that smaller droplets 
defocus faster than larger ones for equal distances from the 
imaging optics plane of focus. Therefor to account for a 
varying sample volume, the droplet size counts for the 
sample volume is multiplied by the ratio of depth of fields. 
The ratio is defined as the maximum measured size's depth of 
field divided by each size bin's depth of field. This 
increases count values for smaller droplet size bins 
exponentially. 

A comparison of sizing the droplets manually vs. 
automatically is given for both experiments in Figures 11 
and 12. Both figures compare the percentage of total count 
vs. droplet size. Also included is the linear, volume, and 
sauter mean diameters for each method. 

Results 

The first experiment is more controlled than the second 
because no decision is necessary on the observer's part as 
to whether or not a sphere (droplet) is in focus, and thus 
all objects were counted except for droplet clusters. In 
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experiment 1, the manual measurement using the video monitor 
varied from the measurement using the micrometer by a 2-6% 
range for the distribution's mean diameters. This includes 
errors induced in the measurement process (selecting the 
object boundaries in both cases is subjective) and actual 


distribution differences 

because of 

the small 

sample size 

(500 for the micrometer. 

175 for 

the monitor) 

and human 

error . 




The distribution's mean 

diameters 

calculated 

using the 

automatic system varied 

over a 

3-9% range 

from the 


micrometer values. It measured 172 objects, 1 of which was 
an error. 

In experiment 2, the manual #2 measurements varied 2-13% 
from manual #1 over the various mean diameters. Because 
this second experiment involved real spray images, a 
decision as to the selection of droplets by their degree of 
focus had to be made by the observers. 400 images were 
analyzed containing approximately 8000 candidate objects. 
The distribution size totalled 150 and 120 droplets for the 
#1 and #2 manual measurements respectively. Both manual 
measurements were averaged to provide mean diameter values 
to compare against the automatic system. The automatic 
system varied 5-16% from those average values. It selected 
85 focused droplets. 
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The variation between measurements in both experiments is 
due to the small sample size and human measurement 
variations. size[18]. 


CLOSURE 

A system has been developed which automatically measures 
droplet size distributions from video images. This offers a 
method to measure areas of spray distributions were 
conventional devices fail such as at the spray nozzle tip. 
This method may also be used to verify the calibration of 
other instruments. The accuracy of this technique has been 
evaluated by two preliminary experiments and was shown to be 
in the range of 2-16%. It is felt that this is very 
satisfactory range. From experience gained in these 
experiments we feel that a next step is a dual purpose 
experiment involving larger sample sizes. This would 
confirm our experimental results and at the same time 
characterize a spray. 
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Figure 1. Equipment block diagram 
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Original Image After Gausian Filter 


Figure 2. Gaussian filtering operation 
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Figure 3. Original Image 
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Figure 4. Estimation of d(x,y) 
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Figure 5. Normalized image 
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Figure 6. Threshold image 


633 


Out of Focus Droplet 

0O / O 

(g°° 


o 


CP 

•o 


V 


O 


© 


rs. 




<o 

o 


o 



Threshold Anomaly 


As*—' — eo^s o o O 

o 

O 



o® 


Mosaic Edge Anomaly 




O 

<0 


o 


O 


Droplet Cluster 

4 

P O 


o 


A 


Figure 7. Boundary image 
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Figure 8. Radial and Standard deviation curve for candidate 
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Figure 9. Quadrant plots for deviation derivatives 
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DROPLET DIFFRACTION INTENSITY FUNCTION 
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Figure 10. Feature measurement 
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MEASUREMENT COMPARISON 
TV-MANUAL vs. MICROMETER vs. AUTOMATIC 
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Figure 11. Synthesized distribution on microscope slide 


638 





MEASUREMENT COMPARISON 
MANUAL vs. AUTOMATIC 
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Figure 12. Heptane spray 
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Figure Legend 


1. Equipment block diagram 

2 • Gaussian filtering operation 

3. Original Image 

4. Estimation of d(x,y) 
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6. Threshold image 

7. Boundary image 

8. Radial and Standard deviation curve for candidate 

9. Quadrant plots for deviation derivatives 

10. Feature measurement 

11. Synthesized distribution on microscope slide 

12. Heptane spray 
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